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We numerically and theoretically study macroscopic properties of dense, sheared granular mate- 
rials. In this process we first introduce an invariance in Newton's equations, explain how it leads 
to Bagnold's scaling, and discuss how it relates to the dynamics of granular temperature. Next we 
implement numerical simulations of granular materials in two different geometries-simple shear and 
flow down an incline-and show that measurements can be extrapolated from one geometry to the 
other. Then we observe non-affine rearrangements of clusters of grains in response to shear strain 
and show that fundamental observations, which served as a basis for the Shear Transformation Zone 
(STZ) theory of amorphous solids can be reproduced in granular materials. Finally we present 

constitutive equations for granular materials as proposed in y], based on the dynamics of granular 
temperature and STZ theory, and show that they match remarkably well with our numerical data 
from both geometries. 



I. INTRODUCTION 

Historically reserved to the engineering community 0, 
H, HQ, g ranu lar materials recently emerged as a new 
field of study for physicists [E II E& El El El Ei : a 
state of matter that is not classifiable according to the 
traditional ternary-solid, liquid or gas- and requires sci- 
entists to rethink the foundations of statistical physics 
and thermodynamics fljl ITtij . The first theories of gran- 
ular materials were motivated primarily by the need 
to predict the creep motion of soils and their stability 
properties. Inspired by continuum theories of plastic- 
ity, these approaches were restricted to quasi-static de- 
formation and incipient failure HSU E3 EH E3- 
In physics, initial progress came in understanding dilute 
granular materials with the development of kinetic the- 
ory m m m m m m m m m m mmm. But 

kinetic theory is based on the assumption that interac- 
tions between grains occur during instantaneous binary 
collisions, a condition that is realized by very few "rapid" 
flows outside the laboratory. 

Traditional approaches thus leave us with little under- 
standing of the elementary mechanisms of deformation in 
granular materials. It is no wonder that the elaboration 
of physically inspired constitutive equations for granu- 
lar materials is currently facing a collection of controver- 
sial issues which challenge its most fundamental aspects. 
Common points of contention include: (i) the relevant 
features of the grain-grain interaction-e.g. Hertzian ver- 
sus hard-sphere repulsion- 0, IS |3; (ii) the domain 
of applicability of kinetic theories, |29l l3ll Ejjl Ejjl I37I : 
(iii) the observation and micromechanical origin of Bag- 
nold scaling in dense flows [S IS IS III El E3 ; (iv) 
the plausible need to formulate "non-local" constitutive 
equations li l |45 Jiti |i ? . |il due to the existence of 
chain forces |4a~50l l5ll~ 52L l53|. Uncertainties about 
these questions have fueled a wealth of theoretical mod- 



els. Some models posit that frustrated rotation plays 
a preeminent role in jamming [Hil l55j . Many mod- 
els propose extensions of kinetic theory by using either 
granular temperature 1561 . introducing strongly density- 
dependent viscosities |57l l58l IHsf . or incorporating a 
quasi-static stress at the internal friction angle of the 
material [S IS IS IS IS EH • Other models are based 
on the introduction of chain- forces Ii3. liH lit^ . activated 
processes [S], "granular eddies" |47| . coexisting liquid 
and solid microphases in a Ginzburg-Landau formula- 
tion [SISISlj or 'spots' of free- volume associated to 
cooperative diffusion |7lLl7a|. 

Several of the above-mentioned models conjecture 
jamming mechanisms-frustrated rotations, chain forces, 
granular eddies-which are so specific to granular materi- 
als that they do not allow for connection with jamming 
in other materials. Our approach rests on the viewpoint 
that, since jamming is observed in many amorphous sys- 
tems, it is likely that it has a common origin. It was thus 
proposed Q that a rescaling of the dynamics of perfectly 
hard granular materials would make it possible to map 
their properties onto more typical glass formers. Building 
on the analogy between granular materials and metallic 
glasses then provides a local mechanism of jamming on 
the basis of the Shear Transformation Zone (STZ) theory 
of plasticity H@ 

The goal of the present paper is both to discuss some 
of the microscopic assumptions underlying the construc- 
tion of constitutive equations for granular materials and 
test the specific predictions of the STZ theory formula- 
tion of such equations Q . Our approach is based on the 
following expectations concerning the four points of con- 
tention identified above: (i) major properties of dense 
granular flows can be captured by perfectly hard grains; 
(ii) jamming involves variations of the frictional stress 
in regimes where the collisional contributions- those pre- 
dicted by kinetic theory- are negligible; (iii) for flows of 
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perfectly hard grains, Bagnold's scaling is relevant at all 
densities and arises from an invariance in the equations of 
motion; and (iv) the rheology of dense granular materials 
is local in a sense to be defined further. 

Indeed, (i) a recent work by Campbell has shown- 
convincingly to us-that most natural and experimental 
flows occur in a regime of strain rates where grains can be 
considered as hard-spheres [Tij. Although this question 
may still be debated, we will take this for granted and re- 
strict our study to a material composed of perfectly hard 
grains, (ii) The breakdown of kinetic theory in dense 
flows has now been characterized in several numerical 
and experimental studies and in- 

dicates that jamming is associated with changes in the 
contribution of long-lasting contacts to the stress tensor, 
not the collisional contribution. We will take this idea 
for granted in the present work, but we do believe that a 
quantitative analysis of this statement is needed and will 
devote a future paper to this question, (iii) We will show, 
following Q that in the case of perfectly hard grains, 
Bagnold's scaling holds for any density. Assumption (iv) 
is strongly inspired by the works by Aranson, Tsimring 
and coworkers, who have shown that a Ginzburg-Landau 
formulation of a fluid-solid mixture j^l E3 could ac- 
count for important properties of granular flows. This 
suggests that a local formulation of granular rheology, 
coupled to hydrodynamic equations, can be sufficient to 
account for a large part of the phenomena observed in 
dense granular flows. However, no numerical study has 
been specifically devoted to this question. 

Therefore, we will first study whether dense granular 
flows can be described as a local phenomenon, governed 
by local hydrodynamic equations, or whether granular 
flows must be treated non-locally, relying on the emer- 
gence of long-range correlations through mechanisms 
such as force chains. For this purpose we have imple- 
mented Contact Dynamics simulations of sheared gran- 
ular materials in two different configurations: (i) simple 
shear in a periodic cell and (ii) thick flows of granular 
material down an incline plane. We will show that the 
measurements obtained in either configuration can be ex- 
trapolated to the other. 

Next we follow Falk and Langer 0, 0] and test for the 
consistency of the STZ picture of material deformation 
advocated by these authors. Our observation provides 
further grounding for the analogy between granular ma- 
terials and common glass-formers and direct support for 
the relevance of STZ theory to granular materials. A re- 
view of the STZ constitutive equations for granular flows 
will follow and finally a fit of the theory with our numer- 
ical data. 

The organization of our paper is as follows. In Sec- 
tion II we present equations of motion in the mathemat- 
ical limit of perfectly hard grains. We discuss the in- 
variance properties of the equations of motion and how 
these properties are related to Bagnold's scaling. We also 



discuss why the results for perfectly hard grains are ap- 
plicable to natural and experimental granular materials. 
In Section III we construct our numerical test, provide al- 
gorithmic details on the Contact Dynamics method, and 
compare the rheology of a dense granular materials in a 
periodic cell and down an incline plane. Finally, in Sec- 
tion IV, we present Falk and Langer's STZ theory, show 
how it adapts to granular materials, and conclude with 
fits of our data. 



II. FUNDAMENTAL RESULTS FOR PERFECTLY 
HARD GRAINS 

When grains are dry-so that no water bridges induce 
attraction-and of size larger than the micrometer scale- 
so that no electrostatic interaction intervenes-their inter- 
action is purely repulsive. The interaction results from 
the elastic deformation of grains at contact and the dis- 
sipation of energy via friction and collisions. The com- 
plexity of this interaction motivates our first question: 
which properties of the grain-grain interaction contribute 
to any particular macroscopic observation? In some in- 
stances, details of the grain-grain interaction seem criti- 
cal: for example, the Hertzian repulsion 0| is essential 
to understand the acoustic properties of granular materi- 
als [iol l8l| . Numerical implementations of granular ma- 
terials have thus relied on more or less elaborate models 
of the grain-grain interaction [il E^l • 

Here we are concerned with dense flows of granular 
materials and, in particular, flows down inclines as found 
in the experiments by Pouliquen [s3 | and numerics by 
Silbert and coworkers 0, 03- For these dense granu- 
lar flows a recent study by Campbell helps us assess the 
importance of the elastic (soft) part of the repulsive po- 
tential El versus the limit in which grains appear as 
perfectly hard spheres. Campbell presented a detailed 
analysis of the different flow regimes obtained in a three 
dimensional simple shear simulation of dense granular 
flows when varying the stiffness k of the repulsion, the 
shear rate 7, and the mass density (j). He found that 
the dimensionless parameter T = ^ybp-, where D is the 
grain size, dictates the character of the flow. This quan- 
tity is directly related to a Mach number which involves 
the ratio of the shear velocity D 7 over the sound speed 
c s : M = Dj/c s — l/VT. The hard-sphere limit cor- 
responds to the situation where sound waves travel very 
fast compared to the motion induced by the shear flow. 
This is the limit of very small Mach number, or small 
shear rates. Specifically, in the numerics by Campbell, 
this limit is reached for Mach numbers below ~ 10~ 2 . 
Since the sound speed in granular materials is of order 
lOOm/s, if we assume a grain size of order 1mm, the 
Mach number is expressible as M = 10~ 5 7. Therefore, 
in order to be in the limit where grains behave as if they 
are perfectly stiff, it suffices to restrain oneself to shear 
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rates below 1000s -1 . Most experimental and natural sit- 
uations occur at shear rates far below this limiting value 
and we can conclude that most flows of granulates are in 
the limit where the soft part of the repulsion is entirely 
masked by the steric exclusion. To study these flows it is 
sufficient to consider properties of perfectly hard grains. 
In the following sections we explore how the mathemati- 
cal limit of perfectly hard spheres gives insight into fun- 
damental processes which relate to most all experimental 
shear flows. 



Equations of motion and Hard Sphere conditions 

The motion of N spherical grains in a d-dimensional 
granular material is determined by Newton's equations 
for the positions angular orientations 6i, momenta pi 
and angular velocities wc 



dqi_ _ Pi_ dpj _ y> 
dt "m,-' dt ~ ^ lj 
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where F ext represents an external force such as gravity, 
Fij represents a contact force on grain i by grain j, hij 
is the unit normal vector pointing from grain i to grain 
j, Ri is the radius of grain i, and Ii is the moment of 
inertia. 

These equations must be complemented with a pre- 
scription for the contact forces. In the hard sphere limit 
these contact forces are determined self-consistently by 
the conditions that (i) there is no penetration between 
grain-a force is instantaneously created upon contact to 
impede penetration and remains non-zero until the con- 
tact is broken-and (ii) by the friction law which couples 
to rotational degrees of freedom. 

Important properties of granular materials arise di- 
rectly from an invariancc of the equations of mo- 
tion Q |2J). We now spend some time studying these 
properties and assessing their consequences for macro- 
scopic observations, in particular Bagnold's scaling. 



Bagnold's scaling 

The success of kinetic theory came in a large part from 
its ability to account for the scaling between stress a and 
strain rate 7 [a ~ 7 2 ) first observed by Bagnold in dense 



granular materials 38] . In Bagnold's experiment, the 
strain rate was set and the stress measured. Bagnold jus- 
tified this behavior by assuming that the stress measured 
in his experiments resulted solely from binary collisions: 
the frequency of collisions and momentum change per 



collision are each proportional to the shear rate, there- 
fore the stress is proportional to the square of the shear 
rate. 

Recent observations have complicated Bagnold's orig- 
inal assessment. On the one hand, the observations by 
Bagnold have been criticized: they may have arisen from 
a secondary instability of the granular flow in his shear 
cell On the other hand, Bagnold's scaling has been 
directly observed by measuring shear stress and strain 
rate profiles in numerical simulations of granular flows 
down inclines 0, > an d is found to be consistent with 
experimental observations of the average flow rate in the 
same geometry [j^. 

Although kinetic theory predicts an exponent which 
agrees with the power law dependence of Bagnold's scal- 
ing, the derivation relies on the assumption that all in- 
teractions between grains occur through binary collisions. 
This assumption is applicable to dilute flows but is not 
upheld for dense flows, such as those originally studied 
by Bagnold. A new argument must be formulated for the 
dense flow regime. 

We contend that in both the dense and "rapid" flow 
regimes, Bagnold's scaling arises from a very fundamen- 
tal invariance of Newton's equations in the hard sphere 
limit. This invariance does not require any of the assump- 
tions of kinetic theory to hold. Namely, for a granular 
material free from external forces, the time evolution will 
obey equations Q 0) with F ext = 0. If we now rescale 
the contact forces by a scalar value F^ — ► F^j A and si- 
multaneously rescale the time t — » ty/A, then Newton's 
equations are transformed to read: 
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where pf ew = pt/VA and uj^ ew = uii/VA. This form for 
Newton's equations is identical to and with new 
values for the momenta and angular velocities. 

Under the rescaling of contact forces and time, the po- 
sitions and angular orientations remain unchanged, while 
the velocities are changed in accordance with the time 
rescaling. If we were to watch a movie of one granular 
flow where the grains have initial velocities pi, uji and 
watch another movie at half the speed where the ini- 
tial velocities are doubled pi — > 2pi, uji — ► 2u>i, the two 
movies would look exactly the same in the hard sphere 
limit. The difference in the dynamics is that the contact 
forces measured in the second movie would be 4 times 
larger than those in the first. 

This invariance is a property of perfectly hard grains 
which must hold in the inertial regime. This includes 
the regime infinitely close to jamming, where multi-body 
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interactions dominate collisional terms and the basic as- 
sumptions of kinetic theory fail. In an experiment, this 
scaling breaks down only when it is no longer appropri- 
ate to model the experimental system by perfectly hard 
grains. Relying on the arguments of Campbell intro- 
duced earlier [73j , we can conclude that most experimen- 
tal granular flows are in the regime where it is appropriate 
to model the system by hard grains. 

This invariance also holds for any value of the resti- 
tution and friction coefficients. The separation of time 
scales for dissipation and shear is not a necessary con- 
dition for the invariance to be upheld. The reason is 
that dissipation occurs at contacts between grains: dis- 
sipative processes are correlated to the motion of grains, 
hence the rate of dissipative processes scales with 7. 

Granular Temperature 

Because the positions and angular orientations of 
grains are invariant to a simultaneous change in time and 
force scales, the path that a granular material takes in 
configuration space is also invariant: only the speed along 
this path is altered. Therefore, the path that a granular 
material takes in configuration space can be separated 
from the rate at which events occur along that path. 

This observation leads to a natural definition of granu- 
lar temperature T: the rate at which microscopic events 
occur on the configuration space path is defined as 
s/T I (R) (the average grain size (R) is inserted so that T 
has units of velocity squared) . These microscopic events 
are not necessarily collisions, but can also be understood 
as force fluctuations or minute displacement of grains al- 
lowing for propagation of "hard-sphere" noise. 

The invariance is directly related to the structure of 
Newton's equations and the existence of inertial terms- 
or accordingly, to the fact that kinetic energy, no matter 
how small, is well defined. The granular temperature de- 
fined above as a "velocity" along phase space trajectories 
is thus naturally proportional to kinetic energy: 

r, = ^((« 2 )-(«> 2 ) + ^((- 2 >-(-> 2 ), (5) 

where m is the mass, v the velocity, / the moment of 
inertia, u> the angular velocity, and brackets denote an 
average over grains. Since the granular and kinetic tem- 
peratures are defined up to a constant factor, we will 
equate them. This means that at any time, the kinetic 
energy provides an estimate of the frequency of elemen- 
tary events in a multicontact system. 

What quasi-static limit? 

Consider a uniform granular material in simple shear 
flow. If we imagine simultaneously changing the force 



and time scales then the invariance in Newton's equations 
guarantees that the dynamics will remain unaltered, re- 
sulting only in new velocities that are determined from 
the time rescaling. However, under the rescaling, quanti- 
ties such as the pressure and strain rate will take different 
values. The quantities of interest are those that are in- 
variant to the transformation. These quantities must be 
formed of ratios of force scales, ratios of time scales, or 
ratios of force to time scales. 

Denote p as the pressure, a the shear stress, 7 the 
strain rate, <j) the mass density, and D the average diam- 
eter of grains. Some important invariant (and dimension- 
less) quantities are: -j^jp-yz , ^ 253 jp, , ^f— > - t > ^75^2^2 j 

and — . These must be single valued functions of den- 
p 

sity only (independent of the strain rate). This observa- 
tion, established by the invariance in Newton's equations 
which holds for granular flows in the dense and collisional 
regimes, automatically predicts Bagnold's scaling: since 
t £ is a function of density only, it follows that a oc j 2 . 

This leads to a simple yet remarkable property: at a 
given density, changing the strain rate does not lead the 
system toward any different regime or "phase" . The ma- 
terial does not go closer to jamming because the strain 
rate is scaled down. The system is always exploring the 
same trajectories in phase space, but at a slower speed. 
In short, there is no quasi-static limit for perfectly hard 
grains. This is apparent in the jamming phase diagram 
of granular materials |8fij where changing the density al- 
lows the granular material to jam, while changing the 
strain rate does not. This property is also directly re- 
lated to the observation by Campbell that "there is no 
path between inertial (rapid) flow and quasi-static flow 
by varying the shear rate at fixed concentration" [7flj . 
This is no surprise once we understand the invariance of 
Newton's equations for perfectly hard grains: only two 
regimes are accessible, jammed or inertial, and density is 
the only parameter which controls the statistical proper- 
ties of a flow of perfectly hard grains. 

III. CONSTRUCTION OF A NUMERICAL TEST 

In constructing a numerical test our goals are to mea- 
sure stress-strain relations when granular temperature 
and stresses are appropriately scaled, show that they 
compare well with the standard response of yield stress 
liquids, and show that the local rheology measured in 
simple shear flow matches the bulk rheology of a granu- 
lar flow down an incline. 

In order to address these issues, we implement numer- 
ical simulations of granular materials in two different ge- 
ometries: 

• We implement simple shear flow in a cell with Lees- 
Edwards (LE) boundary conditions. In this config- 
uration, the density and the shear rate is prescribed 



FIG. 1: Snapshot of a granular material simulation in the 
simple shear configuration. Each grain has an average velocity 
in the x-direction given by yy, where 7 is the strain rate. The 
center of the cell is defined as x — y = 0. 

and the simulation cell is, by construction, transla- 
tionally invariant. This grants direct access to av- 
eraged quantities of the granular temperature and 
stress tensor. Using this configuration we can char- 
acterize the steady state relation between stresses, 
granular temperature and strain-rate and extract 
numerically the parameters of a constitutive law 
for granular materials. A screenshot of this shear- 
ing geometry is shown in Figure ^ 

• We implement granular flow down an inclined plane 
made of stationary grains. The simulation cell is 
periodic in the direction (x) parallel to the plane 
and the flow is inhomogeneous in the perpendicular 
(y) direction. In this configuration, the stresses are 
prescribed by the angle of the incline. We perform 
x-averaged, y-dependent measurements of granu- 
lar temperature, velocity profiles and strain-rate. 
Large heights of the granular layer grant access to 
the bulk rheology of the flow. This permits us to 
check the existence of a well-define bulk rheology in 
the large height limit, and to compare it with the 
measurements in simple shear. A picture of this 
shearing geometry is shown in Figure 

In order to make a quantitative comparison of the two 
simulations, we use the same material: a two-dimensional 
polydisperse mixture of constant density grains with the 
radii drawn from a flat distribution with average radius 
(R) and width a. For all of the simulations in this paper 
we set a/(R) = 0.5, using (R) = 0.7. This distribu- 
tion prevents crystallization and produces an amorphous 
granular material, as can be seen from measurement of 
the pair correlation function in Figure [3J The horizontal 



FIG. 2: Snapshot of a granular material simulation in the 
incline flow configuration. Fixed grains (indicated by filled 
circles) create a stationary incline at angle 8 on which the 
flowing grains are accumulated and allowed to flow. Gravity 
drives the motion and is directed vertically downward. 

axis (d) in this figure corresponds to the distance between 
a pair of grains, divided by the sum of their radii. This 
normalizes the figure so that d = I corresponds to con- 
tacting grains. Other than this peak, the function has 
some small variation (which implies a correlation) be- 
tween d = 1 and d = 3. However, there is no correlation 
beyond d — 3. Because there is no large scale correlation, 
this implies that the granular material is amorphous. 

In this paper we present data for normal and tangential 
coefficients of restitution given by e n = et = 0. It was 
shown by Chevoir et al 86] that the regimes reached by 
granular materials are almost independent of the restitu- 
tion coefficients below some threshold around e n = 0.7. 
Different friction coefficient have been used: frictionless 
grains and /j, = 0.4. 

Our simulations of granular materials rely on the Con- 
tact Dynamics algorithm [13, EM E3, H3, lEll We refer 
to the literature for technical details and present below a 
brief overview of the method, supplemented by a single 
addition we made in order to construct a Lees-Edwards 
simulation cell for frictional granular materials. 



Contact Dynamics 

A Contact Dynamics algorithm was constructed to 
carry out numerical simulation of spheres interacting 
through the enforcement of hard sphere conditions [§3 • 

When a contact occurs there is an non-continuous force 
created that prevents the contacting grains from pene- 
trating. The magnitude of this force is chosen to ensure 
that the final relative velocity u of the grains is related 
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FIG. 3: The pair correlation function for a frictionless gran- 
ular material at high density in linear (top) and log (bottom) 
scales. This function is representative of the pair correlation 
functions for other densities and for granular materials with 
friction between grains. There is a large peak correspond- 
ing to contacting grains (d = 1) and some variation between 
d = 1 and d = 3. For d > 3 there is no correlation between 
grains. This implies that the material is amorphous. 

to the initial relative velocity u° via the equations 

u„ = -e n u° n , u t = e t u°, (6) 

where e n , et are constant restitution coefficients that will 
depend on the shape and consistency of the grains, and 
the n, t subscripts represent the normal and tangential 
parts of the relative velocity with respect to the contact. 
At each time step, the algorithm computes the contact 
forces by first ensuring that these relations hold at each 
contact. To include friction, the Contact Dynamics algo- 
rithm ensures that the resulting tangential force F t is less 
than or equal to [iF n where /i is the friction coefficient 
between grains and F n is the normal force. If this con- 
straint does not hold, then the algorithm sets F t = fiF n 
in order to comply with Coulomb friction. 

In this way the Contact Dynamics algorithm calcu- 
lates, at each time step, contact forces that are consis- 
tent with Newton's equations and the hard sphere con- 
tact law. 



SLLOD equations for simple shear flow 

Lees-Edwards (LE) boundary conditions permit us to 
prescribe the deformation of a material by controlling the 
positions of the image cells • In all of the simulations 
presented here, we impose a constant strain rate 7 so 
that a grain at position y has an average velocity of jy 
in the x-direction (see Figure 0. 



It was recognized in early implementations of LE 
boundary conditions that when deformation is applied 
through the image cells, the information needs time to 
propagate from the cell boundaries to its center. In order 
to ensure rapid propagation of this information and pre- 
vent the boundaries between cells from making unphys- 
ical contributions to the motion, it is necessary to mod- 
ify Newton's equations by introducing so-called SLLOD 
terms. These terms can be understood as a sort of "shear 
bath" , with all particles in the cell being directly coupled 
to the overall deformation In practice, the SLLOD 
terms introduce a mechanical perturbation to the equa- 
tions of motion that gives each grain an average velocity 
consistent with simple shear flow. If we separate the mo- 
mentum pi of each grain i into the average part m^yi 
and fluctuating part pi, so that pi — m^yi + pi, then the 
SLLOD equations read: 

The equation for the position qi is simply the result of 
writing the momentum in terms of an average and fluc- 
tuating part. The equation for pi contains a new term 
xj^pi ■ y) which forces the shear flow. Since every grain in 
the primitive cell is acted upon by this mechanical force, 
the constant strain rate is imposed on all of the grains 
simultaneously at the beginning of the simulation. Fur- 
thermore it can be proven that, in the LE geometry, the 
SLLOD equations give an exact representation of simple 
shear flow arbitrarily far from equilibrium |93t l94j . 

For a granular material with non-zero friction coeffi- 
cient /j,, the equations of motions should incorporate ro- 
tations of the grains (for fi = the tangential contact 
force is always zero and there is no rotation). It is ex- 
pected that a SLLOD term should arise in the equations 
of motion for the angular velocity since, in the linear ve- 
locity profile indicative of simple shear flow, the top and 
bottom of every grain should be moving with slightly dif- 
ferent velocities. This will give each grain an average ro- 
tation of 7/2 which must be incorporated in equation (0) 
just as the average velocity xj(qi • y) was incorporated in 
equation JJJ. This leads to the following equations: 

1 j 

where H>i denotes the fluctuating part of the angular ve- 
locity and we have inserted the moment of inertia of con- 
stant density disks in two dimensions. 

Equations J7J and JSJl now give an exact representation 
of simple shear flow for a frictional granular material ar- 
bitrarily far from equilibrium. 

The primary interest of this procedure is that it per- 
mits us to simulate a sheared granular material with a 
homogeneous shear rate. Experimental procedures, e.g. 
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in a Couette cell, do not guarantee that the strain rate 
is homogeneous: the existence of walls induce a non- 
uniformity of the flow and possibly localization of the 
deformation. Our protocol grants direct access to the 
rhcology of the granular material in a self-averaging sit- 
uation. 



Macroscopic Quantities 

We define the stress tensor T,"^ 3 via Cauchy's equation: 
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where (v a ) is the average velocity in the a direction (aver- 
aged over all grains), <f> is the mass density, and f ext is the 
external force per volume. This equation simply states 
that the total time derivative of the average velocity (left 
hand side) is proportional to the divergence of the stress 
tensor, plus any external forces. Cauchy's equation gives 
a definition of the stress tensor from which there exists a 
procedure to derive the functional form of the stress ten- 
sor. This is called the Irving-Kirkwood derivation [94| 
and it yields 
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(10) 



where Vi is the fluctuating velocity of grain i determined 
by vf = v™ — (vf), and V is the volume of the granular 
material (or area in two dimensions). This symmetric 
stress tensor can be written in terms of three variables: 
the pressure p, shear stress a, and first normal stress 
difference N\ defined as 



(11) 



where the signs are chosen so that shear stress and pres- 
sure are positive in our conventions. 

The granular temperature T is measured as 



T 



V 2 



E^ 



(12) 

which is proportional to the kinetic temperature from 
equation |J5}, with the factor of 1/2 coming from the mo- 
ment of inertia calculated for constant density discs. 

Lastly, for all of the numerical data that will be pre- 
sented, we quantify the density of the system by its pack- 
ing fraction v. The packing fraction is defined as the area 
occupied by grains divided by the total area of the sys- 
tem. In our simulations, packing fraction is proportional 
to both the mass density <j> (tt <fi = Av ) and the number 
density n (nir(R 2 ) =v)- 
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FIG. 4: Raw data of the normalized pressure as a function 
of total shear for two frictionless granular material with pack- 
ing fraction 0.8. Data from simulations with different shear 
rates 7 are shown. The top plot corresponds to 7 = 10~ 2 
and the bottom to 7 = 10 2 . The pressure is normalized by 
•y 2 which collapses the two data sets on to one master curve 
(i.e. with this rescaling the top and bottom traces appear 
essentially identical, up to numerical noise), as predicted by 
the invariance for hard sphere systems. 



IV. TEST OF LOCAL RHEOLOGY 

We first study the average relation between pressure p, 
shear stress a, shear rate 7, and granular temperature T 
in simple shear, using the periodic and translationally in- 
variant LE cell. Next we compare these results with data 
obtained for the granular flow down an inclined plane. 

Simple Shear 

Preliminary test 

In all of the simple shear simulations presented here 
we have simulated 2500 grains in a square primitive cell, 
although we have conducted a limited number of simu- 
lations with up to 10000 grains to ensure the accuracy 
of our observations. Because the contact dynamics al- 
gorithm induces some amount of numerical noise, the 
motion of a collection of grains driven at different shear 
rates is not expected to reproduce exactly the same phase 
space trajectory. In Figure^we show raw data of the nor- 
malized pressure pj~ 2 as function of shear strain (strain 
rate multiplied by time), at a packing fraction of 0.8 with 
no friction and at two different values of the shear rate. 

According to the invariance in Newton's equations 
P7~ 2 should be independent of 7, and this behavior is 
confirmed by the measurements in Figure Although 
the shear rates in the two plots differ by a factor of 10 4 , 
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FIG. 5: Invariant quantities a/p (top), j/VT (middle), and 
T/p (bottom) as a function of shear strain for a frictionless 
granular material at packing fraction 0.8. 



FIG. 6: Steady state values of pj~ 2 (squares), a^~ 2 (dia- 
monds), T7~ 2 (circles), and a/p (stars) as a function of pack- 
ing fraction for a frictionless granular material. 



the normalized pressure is virtually identical for both sys- 
tems. Interestingly, not only do the steady state values 
show no shear rate dependence, but the initial transient 
is virtually identical for both values of 7. The invari- 
ance in Newton's equations also predicts that aj^ 2 and 
Tj~ 2 are independent of 7. For all of the simulations we 
have carried out these predictions from the invariance are 
upheld- although numerical noise often disrupts the per- 
fect invariance for large values of shear strain, we see no 
change in the steady state values of normalized pressure, 
shear stress, or granular temperature as the shear rate is 
varied at constant density. These tests ensure that the 
simulations we conduct respect the invariance in New- 
ton's equations for hard spheres. 

For a granular material characterized by its pressure p, 
shear stress a, temperature T, and strain rate 7, we can 
construct three independent invariant quantities: cr/p, 
7/VT, and T/p. In Figure El we show values of these 
three independent invariant quantities as a function of 
shear strain for a frictionless granular material at packing 
fraction of 0.8. For all quantities, steady flow is reached 
by a shear of approximately 0.5, and we will subsequently 
provide stationary data by time-averaging our measure- 
ments between strains of 2 and 10. The values of a/p 
and T/p fluctuate much more than 7/VT. This is due 
to the fact that a and p depend on the forces between 
grains, which are highly fluctuating in the hard sphere 
limit. In our simulations with 10000 grains we observe 
that the fluctuations decrease while the average value re- 
mains constant. This suggests that in the limit of large 
system size, the fluctuations would disappear. 



Liquid-solid transition 

In Figure^lwe present the steady state values of pj~ 2 , 
(77 ~ 2 , Tj~ 2 , and a/p in simple shear for a range of high 
packing fraction systems that we have studied, at zero 
friction. Although there is relatively little change in these 
quantities for small packing fraction, for packing frac- 
tions larger that 0.75 there is a large increase in the values 
of the stresses and granular temperature. Additionally, 
the functional form of the stresses and granular temper- 
ature changes from approximately exponential to a func- 
tion that grows faster than an exponential at v ss 0.75. 

Alam and Luding [95| have measured the steady state 
values of the first normal stress difference Aq , defined via 
equation l|ll|) , in simple shear flow using a monodisperse 
collection of grains. In dilute flows, a non- vanishing value 
for Ni results from collisional terms and the anisotropy in 
the distribution of velocities at Burnett order (9(| . Alam 
and Luding reported that Aq becomes negative at the 
onset of crystallization. We have thus measured the first 
normal stress difference in our system in order to ensure 
the absence of crystallization and as a signature of the 
breakdown of kinetic effects. The observation in Figure[7| 
that Ni > for all packing fractions in our system is 
consistent with the observation that our system remains 
amorphous. The decay of N\ over all packing fractions is 
consistent with the idea that kinetic effects become less 
important as the packing fraction increases, especially 
after v 0.75 when N\ begins to quickly decay. 

In our numerical simulations, we expect the collisional 
contributions to the stress tensor to be negligible. A de- 
tailed study of the crossover between collisional and non- 
collisional regimes will be the topic of a future work (9?| • 
For now we rely on the observation that the kinetic ef- 
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FIG. 7: First normal stress difference N\ as a function of 
packing fraction for a frictionless granular material. TVi > 
for all packing fractions suggests that the granular materials 
are amorphous, and the decay of N\ as a function of packing 
fraction suggests that non-kinetic effects dominate at high 
packing fraction. 

fects, as probed by the first normal stress difference, de- 
cay close to jamming. We also refer to prior works which 
support the same conclusion 

Flow down an inclined plane 

In the simple shear cell the density and strain rate were 
specified, while stresses and granular temperature were 
measured. We now focus on flow down an inclined plane 
which provides a complementary situation. For incline 
flow, stresses are specified by the choice of an angle of 
inclination and by the gravitational field. Then the pro- 
files of velocity and velocity fluctuations are measured 
and grant access to profiles of strain rate and granular 
temperature. 

We report in Figure [S] the packing fraction, granular 
temperature, and A f/y/T as a function of height for the 
steady flow of a non-frictional granular material at an an- 
gle of 12°, with a total height of approximately 50 grains. 
We have conducted simulations with heights ranging be- 
tween 25 and 100 grain diameters to ensure that our re- 
sults do not depend on the size of the system. 

We observe that, in the bulk central region of the flow, 
the packing fraction profile is uniform, the granular tem- 
perature is linear, and j/yT is constant. These obser- 
vations hold in our simulations for all angles where the 
granular flow reaches a steady state. We only use data 
from these steady state flows in this paper. These ob- 
servations are consistent with previous observations by 
Silbert et al HHH. 

There are four quantities of interest for incline flows- 
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FIG. 8: Profiles of packing fraction (top), granular tem- 
perature (middle), and j/VT (bottom) as a function of the 
height (y) in the pile, measured in grain diameters, for a non- 
frictional granular material at a 12° incline. 

p, a, 7, and T— and these lead to three independent in- 
variant quantities A f/^/T, cr/p, and T/p. We observe in 
our simulations that all of these invariants are constant in 
the bulk of the incline flow. Therefore it is legitimate to 
compare these constant values with the constant values 
obtained from simple shear simulations. In Figure we 
present how the constant values of a/p, T/p, and j/yT 
in the bulk of the flow depend on packing fraction and 
compare with our results from the simple shear cell. The 
fact that data from different shear flows fall on the same 
curves is remarkable and suggests that one theory should 
be able to describe simple shear and bulk incline granular 
flow. 

Interestingly, the data from different flows do not over- 
lap over a large interval of packing fraction: the flow 
down an inclined plane provides values at higher values 
of packing fraction than the simple shear cell. This is due 
to (i) steady flows down an incline plane are more eas- 
ily reached for lower inclinations, hence higher densities; 
and (ii) the simple shear deformation is more difficult to 
integrate numerically at higher densities, because the pe- 
riodic cell induces additional constraints that the Contact 
Dynamics algorithm manages with difficulty. Neverthe- 
less, our use of two different configurations grants access 
to a broad range of ^/VT, a/p, and T/p; and the sets of 
data are consistent with the existence of a unique, local 
relation between them as apparent in Figure E3 



Origin of a Local Rheology 

The excellent agreement between the two sets of data 
challenges the belief that non-local effects arise in dense 
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FIG. 9: The values of a/p (top), -y/VT (middle), and T/p 
(bottom) plotted as a function of packing fraction. Data from 
simple shear flow (circles) and flow down an incline (squares) 
match on the same curves. This suggests that there is a lo- 
cal rheology that is independent of the particular shearing 
geometry. 



granular flow. The data in Figure [5] supports a very con- 
servative opinion: the motion of the grains decorrelates 
beyond some finite length scale, in accord with the fast 
decay of the pair correlation function (see Figure |3J) . 

As a consequence, in the bulk of the flow down an in- 
cline, it is possible to view layers of granular materials 
as effective simple shear cells. Such a layer of granular 
material at height y responds essentially as if it was con- 
fined in a simple shear cell, in the absence of body forces, 
with sustained external stresses a(y) and p(y). Of course 
the invariance in Newton's equations, which holds exactly 
for the simple shear cell, is slightly broken by the grav- 
itational force field. However, deep in the bulk of the 
flow, large confining stresses eventually dominate over 
the gravitational field. This approximate invariance suf- 
fices to predict that Bagnold's scaling must hold for the 
bulk regions of incline flows and explains the numerical 
data of Silbert and coworkers [4lt |42| . 

Hydrodynamic equations are expected to arise when 
the gradients of macroscopic variables become small com- 
pared to the macroscopic variables themselves. What is 
surprising, however, is that the locality of the rheology is 
observable for the moderate heights that we can access 
in our numerical simulations. We believe that the reason 
why the hydrodynamic regime seems to be observed at 
accessible scales is the following: (i) the relevance of hy- 
drodynamics limit depends on how a macroscopic quan- 
tity 'i? compares with fV'F, where £ is a length-scale; (ii) 
I is not the size of the grain or a cluster of grains- it is the 
size of the region sampled by a given grain per strain unit, 
which is on the order of the mean-free path. Since for a 
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FIG. 10: a/p plotted against ■y/VT in steady state simple 
shear flow (circles) and incline flow (squares). Filled sym- 
bols correspond to frictional granular materials with coeffi- 
cient of friction y = 0.4 and open symbols correspond to 
non-frictional granular materials with y = 0. 



dense granular material £/(R) « 1, the validity of the 
hydrodynamic limit should hold over moderate heights. 

As we will see, our theory provides constitutive equa- 
tions which relate the static friction coefficient a/p to 
the ratio 7/vT. Anticipating the following sections, we 
present in Figure HOI a plot of a/p versus j/VT, with 
data from both the simple shear and incline flow geome- 
tries. As we will argue, such a plot is expected to be the 
granular counterpart of a stress versus strain rate plot 
for a glassy material Q, with the shear stress normal- 
ized by the pressure and the strain rate normalized by 
the granular temperature. We sec here that the analogy 
is striking: when rescaled properly the granular material 
presents typical features of normal yield stress fluids [§3 • 
For large values of normalized strain rate, the normalized 
stress is proportional to the normalized strain rate. For 
small values of the normalized strain rate, the linear rela- 
tionship no longer holds and there is a yield (normalized) 
stress at zero (normalized) strain rate. 



V. CONSTITUTIVE EQUATIONS 

Several interesting results stem from the preceding ob- 
servations: (i) for a given density, only invariant quanti- 
ties are relevant to describe the state of a granular flow; 
(ii) when invariant quantities are considered, the granular 
material displays a "normal" rheology of a yield stress liq- 
uid; and (iii) the rheology measured in the Lees-Edward 
cell extrapolates to the rheology measured from bulk data 
of incline flow. 

These observations foster our hope to construct local 
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constitutive equations for dense granular materials. Be- 
fore addressing the specificities of the STZ formulation 
of constitutive equations, we wish to argue that the scal- 
ing invariance of Newton's equations calls for a parallel 
scaling form of the constitutive equations. 



Scaling form for constitutive equations 

In general, constitutive equations should provide a re- 
lation between stresses and strain-rate, and a set of state 
variables {ipi} which characterize the internal structure 
of the granular material. For granular materials the gran- 
ular temperature T plays a quite specific role as a state 
variable. Other state variables are expected to account 
only for geometric properties of the granular packing. 



Conservation of energy 

Ogawa was the first to recognize that in granu- 
lar materials the temperature cannot be prescribed by a 
thermal bath, but is set by energy balance. Kinetic the- 
ory provides estimates for the energy dissipation rate in 
dilute systems, hence providing approximations for the 
equation of motion which governs granular temperature. 
As a result, the notion of granular temperature has been 
tied to kinetic theory. 

Granular temperature, however, is not reserved for the 
description of dilute flows. It is relevant in all inertial 
flows since, as discussed previously, there is no quasi- 
static limit for hard-spheres, even infinitely close to jam- 
ming. In particular the dense flows of our simulations 
do not verify the assumptions of kinetic theory, yet in 
these inertial regimes kinetic energy is a perfectly rele- 
vant physical observable and is set by energy balance. 

At any time, the variation of kinetic energy results 
from the balance between external work done on the sys- 
tem and dissipative mechanisms. Because the system ex- 
plores phase-space trajectories at a velocity proportional 
to VT, the rate of energy dissipation should scale as 
•v/TT. The square root gives the frequency of dissipative 
events and each event dissipates an energy proportional 
to T. Energy is introduced to the system via the external 
forcing with a rate of 177. Therefore, the energy balance 
equation should be of the form 

<j>f = aj -a(<f>, {i>i})Tj^, (13) 

with <fi being the mass density. The factor a is a ge- 
ometric factor which depends on the relative positions 
of grains, but should not incorporate any further depen- 
dence on T or the stress tensor. In particular, it should 
depend on the mass density cp and possibly on other state 
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FIG. 11: a from equation 113H as a function of shear strain 
for a simple shear flow of a non-frictional granular material 
at packing fraction 0.8. For small values of shear a varies 
slightly, but quickly becomes constant at a shear strain of 
approximately 0.04. 

variables {ipi} which characterize the geometrical struc- 
ture of the granular packing. By scaling invariance, the 
steady state value of a should depend of cp only. 

This equation has a similar form to the energ y conser- 
vation equation derived in kinetic theory [lOfJ. This is 
no surprise since kinetic theory must uphold the invari- 
ance of Newton's equations. Kinetic theory provides an 
estimate for a(<p) and it is possible that even in dense 
flows this estimate remains reasonable for non-frictional 
grains. However, multibody collisions may induce de- 
parture of this relation from the predictions of kinetic 
theory. 

In Figure ^2 we present a plot of a as a function of 
shear strain for a non-frictional granular material in sim- 
ple shear flow, calculated via equation lfP3*)l . We observe a 
dynamic in a at the start of the simulation which quickly 
disappears. The dynamic is likely due to the effects of 
an unknown state variable ipi. 



Constitutive relation 

Since granular temperature is defined as the frequency 
of elementary events along phase-space trajectories, the 
strain rate should scale as 

^ = /(-,-,{>,}), (14) 
VT \P P J 

where / denotes an unknown function and, once again, 
the state variables {ipi} are purely geometric. 

If we combine equation 11411 with equation (|13|l in 
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steady state where T = 0, then we see that 

-3" 



a / a , 

(R) \p 



7 



(15) 



where F is equal to the function /, with T/p evalu- 
ated using equation l|13fl . This equation implies that if 
constitutive relations are written in the form of equa- 
tions i|13[l and 1)140. then Bagnold's scaling is upheld. 

In the rest of this section we study the STZ formulation 
of constitutive equations and apply it to dense granular 
materials. We will see that the STZ theory, when applied 
to granular materials, makes a prediction for the function 
/ in equation 114(1. 



STZ theory 

Basics 

ie (STZ) Theory of 
amo rphous solids was proposed in Q, 0, llOll ll02L Il03l 
Il04| to account for the behavior of dense amorphous ma- 
terials at low temperature. Th e th e ory i s motivat e d by 
observations fro m sim ulations |l05l Il06l Il07l llOSt Il09| 
and experiments |l!0| which suggest that plastic defor- 
mation in amorphous materials results from non-affinc 
rearrangements of small clusters of particles 



111| . Ad- 



ditionally, Falk and Langer were able to show that there 
exist different types of zones which present a preferential 
response to different orientations of shear forces. They 
introduced the densities of these zones as state variables 
to characterize the internal structure of the molecular 
packing. 

Central to the theory is the assumption that once an 
STZ undergoes an elementary rearrangement in a given 
direction it is unlikely that it can shear again in the 
same direction, although it can easily shear in the re- 
verse direction. Thus zones appear as two state systems, 
the states corresponding to the zone orientation being 
aligned (denoted — ) or anti- aligned (denoted +) with the 
shear stress. A rearrangement corresponds to a transition 
of a zone from a ±-state into a =F-state and vice-versa. 
The plastic shear rate is given by the rate at which STZs 
respond to external stresses: 



7 oc i?+n + — i?_n_ 



(16) 



where R± are the stress-dependent probabilities that 
zones of ± types are transformed into one another. 

This constitutive relation must be complemented with 
an equation of motion for the densities n± , which is pos- 
tulated to be of the form: 



n± = -Rzpftzp — R±n± + w(a — bn±). 



(17) 



The first two terms correspond to the transformation of 
STZs into their two possible states. The last term ac- 
counts for the fact that STZs are renewed by the overall 



macroscopic deformation: it contains a creation and de- 
struction rate, both proportional to the plastic work w 
of external forces per time unit. 



Observation of directional response 

Before applying STZ theory to granular materials, we 
check that the same qualitative observations as in Q can 
be performed in these systems. Namely that non-affinc 
motion occurs in localized regions and that the positions 
of the localized regions depends sensitively on the orien- 
tation of the shear. The first observation motivates the 
choice of density of STZs as a state variable, and the sec- 
ond observation shows that each STZ has an orientation 
and therefore only responds to a certain orientation of 
shear stress. 

These observations are based on a measure of the non- 
affinity of the deformation of a cluster of a few molecules 
or grains. Following [1|, the grains undergoing non-affinc 
rearrangement can be determined by calculating, for each 
grain, the local strain rate at time t — At. Then, by 
measuring the difference D between the actual position 
at a later time t and the position predicted from the 
local strain rate at time t — At, we can determine which 
grains have moved non-affinely in the time period At. In 
practice, D is determined by minimizing 

b\t, At) = J2 £ (<W ~ r o W - £(<% + iij) 

n i j 

x[ri(t-At)-rl(t-At)]y (18) 

with respect to the shear rate 7^, where the indices i 
and j are spatial coordinates and the index n runs over 
all grains within two diameters of the reference grain, la- 
beled by the index n = 0. The minimum value of D, 
denoted D, is an approximation of the local deviation 
from affine displacement for the reference grain in the 
time interval [t — At, t] . If there is no non-affine motion, 
then the motion of each individual grain should be com- 
pletely determined by a local shear rate and D = 0. If 
there is non-affine motion then D > 0. 

We have applied this test for non-affine motion to gran- 
ular materials in simple shear to produce Figur elT^l This 
figure is the counterpart of Figure 3 in 101] and Fig- 
ure 7 in 0, which were created from simulations of an 
amorphous Lennard-Jones solid. Each picture has been 
created by shearing an identical initial arrangement of 
particles in a certain direction. D is the local measure of 
non-affinity obtained by comparison between the initial 
and final states. If D is larger than a reference value, the 
particle is said to have moved non-affinely and colored 
black. 

In Figure H2*l part (a) the system is sheared from strains 
of to 0.05, and in (b) the system is sheared to from 
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FIG. 12: Screenshots of granular materials in steady state 
simple shear flow with grains undergoing non-afhne displace- 
ment colored black. In (a) the system is sheared from to 
0.05, in (b) the system is sheared from .1 to .15, in (c) the sys- 
tem is sheared from to .15, and in (d) the system is sheared 
in the opposite direction, from to —.05. 



strains of 0.1 to 0.15. We notice that there is a tendency 
for the regions of non-affinc displacement to form clus- 
ters, and the size of the resulting non-affine regions is 
about the same in both (a) and (b). In (c) the system 
is sheared from strains of to 0.15 and now the size of 
the non-affine regions increases, suggesting many more 
fundamental rearrangements of STZs in the larger time 
period. In (d) the system is sheared from strains of to 
—0.05 (in the opposite direction), starting from the same 
initial configuration as in (a) . Once again we observe that 
the non-affine regions tend to form clusters. However, in 
comparing (a) and (d), we notice that the size of the non- 
affine regions is about the same in the two figures, but 
the locations are different. If the regions undergoing non- 
affine displacement did not have an orientation we would 
expect non-affine motion to occur in the same location, 
regardless of the orientation of the stress. However this 
is not the case and the data in Figure IT21 suggests that 
the regions that move non-affinely have an orientation. 
This is qualitative evidence that the core assumptions of 
STZ theory are upheld in granular materials. 



The STZ densities n± account for structural properties 
of a molecular or granular packing. They are thus ex- 
pected to depend on the positions of grains, orientations 
and distribution of forces, orientations of velocities, but 
not on the overall amplitude of the forces or amplitude of 
velocities. Following |3j, these functions are determined 
using the invariance in Newton's equations. Since w rep- 
resents the plastic work done on the system per unit time 
it should be proportional to cry. In order to make w in- 
variant, we normalize by pressure so that w — o~j/p. As 
for R±, because of the invariance in Newton's equations 
we can separate the rate at which an STZ attempts to 
rearrange from the probability that an attempt leads to 
a successful rearrangement. The attempt rate must be 
proportional to \ff which sets the microscopic event rate 
and the probability to rearrange is written as an exponen- 
tial activation factor of the invariant form e ±K<T / p . This 
yields R± oc VTe ±Ka ^. 

Combining the expressions for R± and w with equa- 
tions Ijlfijl and 1)17(1 . while making a change of variables 
from n± to A a n_ — n + and A oc n_ + n+, yields the 
following STZ equations for granular materials: 



7 oc v T (A sinh (na/p) — A cosh {ko~ /p)) 



A oc 7 I 1 -(-A 

A oc 7-(l - A). 
P 



(19) 



A denotes the total density of zones and A measures 
the mismatch between zones of different orientations and 
therefore is related to the anisotropy of the granular pack- 
ing, k and £ are constants that do not depend on the 
macroscopic variables 7, T, a, or p. However we would 
expect these constants to depend on properties of the 
grains such as shape, distributions of radii, restitution 
coefficients, friction coefficients, or other local static vari- 
ables, including density. 



Steady states 

These equations present two types of steady state so- 
lutions Q 0- One branch of solutions represents 
a jammed state 7 = 0, and occurs when A/A = 
tanh (kct/p). The other branch of solution represents the 
steady flow and occurs when A = 1 and A = p/(ja). 

An elementary analysis of the phase diagram of this dy- 
namical system indicates that the jammed state is stable 
if and only if 



C- tanh ( k— ) < 
P V Pj 



(20) 
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and the flowing state is stable otherwise. The limit of 
stability occurs at a critical angle 9* , which is the solution 
of 

Ctanfl* tanh(Ktan0*) = 0. (21) 

6* is identified as the repose angle of our granular mate- 
rial. 

In the steady flowing regime 7 > 0, STZ theory yields 
the following relation: 

'if v \ 

— j= oc I sinh(/€<T/p) - cosh (kct/p) I (22) 

In particular, in the limit when the ratio vanishes, 
a/p converges towards tan 9*. Therefore the STZ theory 
accommodates cases where there is a residual pressure 
and shear stress at zero shear rate, and predicts that 
in this case the shear stress will be proportional to the 
pressure. 
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FIG. 13: Density dependence of a, denned through equa- 
tion 113(1 . for frictional (stars) and non-frictional (circles) sim- 
ple shear flows. 



Numerical tests 

The constitutive equations for granular materials in- 
troduced in this paper are entirely specified by equa- 
tions (fTTTft and (fTT)|) . We now compare the steady state 
relations predicted by these equations with our numeri- 
cal data, which provides independent access to a, p, T, 
and 7 in different shearing geometries. 

We first test equation (|13fl in simple shear flow. We 
find that for all densities investigated, much like Fig- 
ure El a is constant as a function of shear strain. Ad- 
ditionally we find that the steady state value only de- 
pends on the density and friction coefficient. We present 
the measured steady state values of a in Figure ED a s a 
function of packing fraction, for frictionless (fj, = 0) and 
frictional (/1 = 0.4) granular materials. For frictionless 
granular materials a appears to vary exponentially as a 
function of packing fraction, whereas for frictional gran- 
ular materials a does not take on an obvious functional 
form. 

Next, we test the steady state STZ prediction from 
equation (|22(l . In Figure El wc have plotted numerical 
data of a/p as a function of j/VT for frictional and non- 
frictional granular materials in both simple shear flow 
and incline flow. The line drawn through the data is a 
fit to equation l|22(l . The fit matches the data from both 
flowing geometries very well . 

To construct the fit, we have used the data from the 
simple shear cell only. This fit permits us to extrapolate 
the rheology for the flow down an incline plane, even 
at the approach of the repose angle. This shows that 
an accurate determination of the coefficients from the 
STZ equation (|22|l in one shearing geometry allows for 
prediction in a different shearing geometry. 
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FIG. 14: a/p plotted against j/VT for frictional (filled sym- 
bols) and non-frictional (open symbols) granular materials. 
The circles correspond to data from simple shear flow and 
the squares to data from steady state incline flow. The line 
is a best fit to the steady state STZ prediction from equa- 
tion 1221 . where only the data from the simple shear flow 
geometry was used to construct the best fit. This shows that 
STZ theory properly predicts the outcome of incline flow ex- 
periments once the parameters (which depend only on grain 
properties) are determined from simple shear experiments. 

VI. CONCLUSION 

We have implemented numerical simulations of dense 
granular flows in order to clarify the microscopic origin 
of jamming and the specific assumptions needed to con- 
struct the STZ formulation of constitutive equations for 
dense granular flows. 
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Our first goal was to obtain detailed information on 
dense granular flows, in particular on the locality of the 
constitutive relation. For this purpose we implemented 
two different simulation configurations: a simple shear 
cell with Lees-Edwards periodic boundary conditions and 
SLLOD equations of motion, and the flow of a granular 
material down a rough incline. In both configurations we 
measured the components of the stress tensor, packing 
fraction, strain rate, and granular temperature. 

We have also observed the normal stress difference and 
shown that it vanishes in dense flows so that, in this 
regime, shear stress and pressure completely character- 
ize the stress tensor. In each configuration we studied, 
the rheology of dense granular materials is entirely speci- 
fied by five quantities: shear stress a, pressure p, packing 
fraction strain rate 7 and granular temperature T. In 
the shear cell v and 7 are prescribed, and in the flow 
down an incline a and p are prescribed. The rheology 
of the granular material requires us to specify three re- 
lations which determine the remaining unknowns. There 
are many equivalent ways these relations can be formu- 
lated, but we wished to formulate them in a way that 
emphasized the invariance of Newton's equations and the 
expected scaling form of constitutive equations. 

We found, as anticipated in that the relation be- 
tween quantities 7/VT and a/p is remarkably similar to 
the stress-strain rate relation in a normal elasto-plastic 
transition. We were surprised to observe the quality of 
the fit to the data of Figure [21 with expression 
without any density dependence of the parameters of the 
theory. We would have expected that these parameters 
might involve strong density dependence. This observa- 
tion strengthens the hypothesis that the distinguishing 
feature of granular materials is the way energy balance is 
prescribed and otherwise, once the variables are appro- 
priately rescaled, they behave like "normal" materials. 

Another relation emerges naturally from the expec- 
tation that energy balance in granular materials should 
take the form in equation (|13|) . The density-dependent 
parameter a accounts for the rate of dissipation of en- 
ergy and we would expect it to be reasonably captured 
by collisional integrals of the type used in kinetic theory 
(at least for non-frictional materials). A proper deriva- 
tion of energy dissipation in dense flows now appears as 
an accessible yet challenging issue that we wish to exam- 
ine in future works. 

Finally, a third relation remains to be specified. This 
relation can be, for example, a prescription for a/p or 
7/vT as a function of packing fraction. We observed 
this relation in Figure but did not attempt to model 
it. This missing relation should arise from some under- 
standing of how the rheology depends on the local den- 
sity. Such a property might arise from studies of free- 
volume formulati ons of the rheology of dense granular 
materials, [7J El EH EH El El] We did not wish to 
elaborate this aspect in order to focus on the stress-strain 



rate relation, which is the most immediate prediction of 
STZ theory. 

As we see, the STZ theory correctly predicts one 
among three relations which are necessary for a com- 
plete description of dense granular flows. This relation, 
coupled to Bagnold's scaling, suffices to account for the 
shape of the velocity profile down an incline. However it 
does not account for the overall amplitude of the velocity 
profile or the uniform value of the density as a function 
of tan0. 

Note that we have compared here data from a Lees- 
Edward cell with the rheology of a flow down an inclined 
plane in the limit where the height of the granular ma- 
terial is large. The agreement between the two sets of 
data indicates that there is no long-range structure which 
governs the flow. This does not mean, however, that no 
gradient or diffusive terms should ever appear in the final 
form of the rheology of granular flows. In particular, we 
expect that a diffusive term should arise in equation 113|l , 
but in the large height limit this diffusive term does not 
appear in the bulk rheology. 

This brings us to the observations by Pouliquen [8J| of 
a relation h s t op {9) which governs jamming. These obser- 
vations have led to the idea that non-local constitutive 
equations might be required to capture the rheology of 
granular materials. However, we expect these observa- 
tions to be "normal" finite size effects. For example, it 
was observed in [3( that an exponentially decaying kernel 
in place of w in equation (|17(l was sufficient to account 
for large values of h stop (8) when — > 0*. The analyti- 
cal form of the kernel used in this work was consistent 
with an exponential decay of correlation in the system, 
again consistent with the notion that short-range cor- 
relations suffice to understand these effects. But other 
theories also successfully introduce finite length scales 
to account for macrosco pic p henomena. For example, 
Bazant's 'spot' model |7l|,[72( relates such a length scale 
to the size of diffusive regions of high free-volume. Such 
concepts are very interesting as they provide routes to un- 
derstanding non-trivial phenomena in granular physics, 
relying on very plausible forms of finite size corrections 
to hydrodynamics. 

Another puzzling observation by Pouliquen was that 
the amplitude of the velocity profile seemed to scale as 
l/ftstop(0)- As we saw previously, this amplitude is gov- 
erned by a(u) and a missing relation between e.g. a/p 
and v. We thus cannot provide an interpretation for these 
observations. 

Our observations on a model system strongly suggest 
that the definition of a local rheology for granular mate- 
rials is, in principle, possible. Moreover, we see emerging 
from our analysis some state variables and their equa- 
tions of motion. This now opens a very exciting route 
toward a set of local hydrodynamic equations for dense 
granular materials: such equations will offer a predic- 
tive tool to further address fundamental and practical 
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questions in numerous situations where flow equations 
for granular materials are much needed. 
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